Chapter 7 Beta diversity

load("data/data.Rdata")
load("data/beta.Rdata")
load("data/data_filtdamr_80.Rdata")
load("data/beta_filtdamr_80.Rdata")
sample_metadata <- sample_metadata_filtdamr
genome_counts_filt <- genome_counts_filtdamr
genome_metadata <- genome_metadata_filtdamr
genome_tree <- genome_tree_filtdamr
genome_gifts <- genome_gifts_filtdamr
phylum_colors <- phylum_colors_filtdamr
beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt_func <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),]
genome_counts_filt_func <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),]
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt_func %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

7.1 Permanova analysis

7.1.1 Strata

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.03602 0.0180094 2.9837    999  0.062 .
Residuals 72 0.43459 0.0060359                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           natural protected urban
natural             0.714000 0.047
protected 0.724845           0.023
urban     0.044422  0.021094      
adonis2(beta_q0n$S ~ type*season, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(region)) %>%
        broom::tidy() %>%
        tt()
tinytable_aqambstjb94pzal7v646
term df SumOfSqs R2 statistic p.value
type 2 2.4713937 0.09990627 4.204275 0.001
season 1 1.0483711 0.04238048 3.566927 0.001
type:season 2 0.9372685 0.03788915 1.594459 0.007
Residual 69 20.2800898 0.81982411 NA NA
Total 74 24.7371231 1.00000000 NA NA
pairwise.adonis2(beta_q0n$S ~ type*season, data=sample_metadata_row, perm = 999, 
        strata = "region")
$parent_call
[1] "beta_q0n$S ~ type * season , strata = region , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.4635 0.08006 4.7898  0.001 ***
season       1   0.9775 0.05347 3.1991  0.001 ***
type:season  1   0.5610 0.03069 1.8360  0.007 ** 
Residual    50  15.2773 0.83577                  
Total       53  18.2792 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.2906 0.07930 4.3514  0.001 ***
season       1   0.6345 0.03899 2.1392  0.002 ** 
type:season  1   0.4088 0.02512 1.3785  0.035 *  
Residual    47  13.9396 0.85658                  
Total       50  16.2735 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   0.9002 0.06603 3.2536  0.001 ***
season       1   0.9641 0.07072 3.4846  0.001 ***
type:season  1   0.4255 0.03121 1.5378  0.040 *  
Residual    41  11.3433 0.83205                  
Total       44  13.6329 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

Neutral

betadisper(beta_q1n$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     2 0.09096 0.045478 5.7232    999  0.006 **
Residuals 72 0.57213 0.007946                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            natural protected urban
natural             0.0120000 0.006
protected 0.0102235           0.772
urban     0.0046867 0.7844374      
adonis2(beta_q1n$S ~ type*season, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(region)) %>%
        broom::tidy() %>%
        tt()
tinytable_iio38vvm9y2joermqdqd
term df SumOfSqs R2 statistic p.value
type 2 2.8331097 0.12162688 5.312689 0.001
season 1 1.1483346 0.04929861 4.306748 0.001
type:season 2 0.9141128 0.03924334 1.714158 0.009
Residual 69 18.3978927 0.78983117 NA NA
Total 74 23.2934498 1.00000000 NA NA
pairwise.adonis2(beta_q1n$S ~ type*season, data=sample_metadata_row, perm = 999, 
        strata = "region")
$parent_call
[1] "beta_q1n$S ~ type * season , strata = region , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.8601 0.10743 6.7297  0.001 ***
season       1   1.0323 0.05962 3.7349  0.001 ***
type:season  1   0.6017 0.03475 2.1770  0.003 ** 
Residual    50  13.8200 0.79819                  
Total       53  17.3142 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.4101 0.08840 4.9486  0.001 ***
season       1   0.7465 0.04680 2.6197  0.001 ***
type:season  1   0.4017 0.02518 1.4097  0.037 *  
Residual    47  13.3922 0.83962                  
Total       50  15.9504 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   0.8917 0.07545 3.8147  0.001 ***
season       1   0.9941 0.08412 4.2528  0.001 ***
type:season  1   0.3486 0.02950 1.4915  0.106    
Residual    41   9.5836 0.81093                  
Total       44  11.8179 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

Phylogenetic

betadisper(beta_q1p$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
Groups     2 0.05728 0.028642 2.661    999  0.067 .
Residuals 72 0.77498 0.010764                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           natural protected urban
natural             0.046000 0.115
protected 0.059462           0.949
urban     0.111009  0.951525      
adonis2(beta_q1p$S ~ type*season, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))) %>% pull(region)) %>%
        broom::tidy() %>%
        tt()
tinytable_knw2hssvrh7i8f5ulrmh
term df SumOfSqs R2 statistic p.value
type 2 0.3561228 0.09721296 4.088057 0.014
season 1 0.1880028 0.05132022 4.316297 0.006
type:season 2 0.1138033 0.03106555 1.306387 0.243
Residual 69 3.0053980 0.82040127 NA NA
Total 74 3.6633269 1.00000000 NA NA
pairwise.adonis2(beta_q1p$S ~ type*season, data=sample_metadata_row, perm = 999, 
        strata = "region")
$parent_call
[1] "beta_q1p$S ~ type * season , strata = region , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2      F Pr(>F)   
type         1  0.29577 0.10235 6.3336  0.016 * 
season       1  0.19224 0.06652 4.1165  0.006 **
type:season  1  0.06695 0.02317 1.4337  0.219   
Residual    50  2.33495 0.80797                 
Total       53  2.88992 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)
type         1  0.12582 0.04567 2.3870  0.225
season       1  0.09019 0.03274 1.7110  0.151
type:season  1  0.06174 0.02241 1.1713  0.367
Residual    47  2.47750 0.89919              
Total       50  2.75526 1.00000              

$protected_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1  0.09601 0.06457 3.2849  0.001 ***
season       1  0.15237 0.10248 5.2131  0.001 ***
type:season  1  0.04012 0.02698 1.3727  0.255    
Residual    41  1.19834 0.80597                  
Total       44  1.48684 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

Functional

betadisper(beta_q1f$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008931 0.0044656 1.4188    999  0.246
Residuals 72 0.226611 0.0031474                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          natural protected urban
natural             0.38300 0.117
protected 0.36301           0.395
urban     0.12598   0.40392      
adonis2(beta_q1f$S ~ type,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))) %>% pull(type)) %>%
        broom::tidy() %>%
        tt()
tinytable_u755npa0j9224kuzne36
term df SumOfSqs R2 statistic p.value
type 2 0.02226324 0.05258339 1.998067 1
Residual 72 0.40112594 0.94741661 NA NA
Total 74 0.42338918 1.00000000 NA NA
pairwise.adonis2(beta_q1f$S ~ type*season, data=sample_metadata_row, perm = 999, 
        strata = "region")
$parent_call
[1] "beta_q1f$S ~ type * season , strata = region , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2       F Pr(>F)    
type         1  0.01858 0.04818  4.2067  0.001 ***
season       1  0.13428 0.34814 30.3994  0.001 ***
type:season  1  0.01199 0.03108  2.7135  0.262    
Residual    50  0.22087 0.57261                   
Total       53  0.38572 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)   
type         1 0.014221 0.04503  3.691  0.002 **
season       1 0.073060 0.23133 18.962  0.004 **
type:season  1 0.047453 0.15025 12.316  0.007 **
Residual    47 0.181089 0.57339                 
Total       50 0.315824 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df  SumOfSqs       R2       F Pr(>F)
type         1 -0.001311 -0.00992 -0.4856  0.863
season       1  0.014945  0.11305  5.5365  0.121
type:season  1  0.007891  0.05969  2.9233  0.192
Residual    41  0.110674  0.83717               
Total       44  0.132200  1.00000               

attr(,"class")
[1] "pwadstrata" "list"      

7.1.2 Without strata

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.03602 0.0180094 2.9837    999  0.063 .
Residuals 72 0.43459 0.0060359                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           natural protected urban
natural             0.721000 0.036
protected 0.724845           0.028
urban     0.044422  0.021094      
betadisper(beta_q0n$S, sample_metadata$season) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.045726 0.045726 15.566    999  0.001 ***
Residuals 73 0.214440 0.002938                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           autumn spring
autumn             0.001
spring 0.00018133       
adonis2(beta_q0n$S ~ type*season, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_bui1gj85apmrl9jc3c61
term df SumOfSqs R2 statistic p.value
type 2 2.4713937 0.09990627 4.204275 0.001
season 1 1.0483711 0.04238048 3.566927 0.001
type:season 2 0.9372685 0.03788915 1.594459 0.009
Residual 69 20.2800898 0.81982411 NA NA
Total 74 24.7371231 1.00000000 NA NA
pairwise.adonis2(beta_q0n$S ~ type*season, data=sample_metadata_row, perm = 999)
$parent_call
[1] "beta_q0n$S ~ type * season , strata = Null , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.4635 0.08006 4.7898  0.001 ***
season       1   0.9775 0.05347 3.1991  0.001 ***
type:season  1   0.5610 0.03069 1.8360  0.017 *  
Residual    50  15.2773 0.83577                  
Total       53  18.2792 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.2906 0.07930 4.3514  0.001 ***
season       1   0.6345 0.03899 2.1392  0.005 ** 
type:season  1   0.4088 0.02512 1.3785  0.100 .  
Residual    47  13.9396 0.85658                  
Total       50  16.2735 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   0.9002 0.06603 3.2536  0.001 ***
season       1   0.9641 0.07072 3.4846  0.001 ***
type:season  1   0.4255 0.03121 1.5378  0.049 *  
Residual    41  11.3433 0.83205                  
Total       44  13.6329 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

Neutral

betadisper(beta_q1n$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     2 0.09096 0.045478 5.7232    999  0.005 **
Residuals 72 0.57213 0.007946                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            natural protected urban
natural             0.0090000 0.007
protected 0.0102235           0.790
urban     0.0046867 0.7844374      
betadisper(beta_q1n$S, sample_metadata$season) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0133 0.0133019 2.1964    999   0.13
Residuals 73 0.4421 0.0060562                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        autumn spring
autumn          0.137
spring 0.14264       
adonis2(beta_q1n$S ~ type*season, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_xyqd38pfj5e2dcczs211
term df SumOfSqs R2 statistic p.value
type 2 2.8331097 0.12162688 5.312689 0.001
season 1 1.1483346 0.04929861 4.306748 0.001
type:season 2 0.9141128 0.03924334 1.714158 0.010
Residual 69 18.3978927 0.78983117 NA NA
Total 74 23.2934498 1.00000000 NA NA
pairwise.adonis2(beta_q1n$S ~ type*season, data=sample_metadata_row, perm = 999)
$parent_call
[1] "beta_q1n$S ~ type * season , strata = Null , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.8601 0.10743 6.7297  0.001 ***
season       1   1.0323 0.05962 3.7349  0.001 ***
type:season  1   0.6017 0.03475 2.1770  0.015 *  
Residual    50  13.8200 0.79819                  
Total       53  17.3142 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   1.4101 0.08840 4.9486  0.001 ***
season       1   0.7465 0.04680 2.6197  0.001 ***
type:season  1   0.4017 0.02518 1.4097  0.116    
Residual    47  13.3922 0.83962                  
Total       50  15.9504 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1   0.8917 0.07545 3.8147  0.002 ** 
season       1   0.9941 0.08412 4.2528  0.001 ***
type:season  1   0.3486 0.02950 1.4915  0.084 .  
Residual    41   9.5836 0.81093                  
Total       44  11.8179 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

Phylogenetic

betadisper(beta_q1p$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
Groups     2 0.05728 0.028642 2.661    999  0.053 .
Residuals 72 0.77498 0.010764                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           natural protected urban
natural             0.039000 0.088
protected 0.059462           0.961
urban     0.111009  0.951525      
betadisper(beta_q1p$S, sample_metadata$season) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01038 0.010385 0.8932    999  0.385
Residuals 73 0.84878 0.011627                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        autumn spring
autumn          0.389
spring 0.34774       
adonis2(beta_q1p$S ~ type*season, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_019avho6222d8y8g0f8i
term df SumOfSqs R2 statistic p.value
type 2 0.3561228 0.09721296 4.088057 0.001
season 1 0.1880028 0.05132022 4.316297 0.003
type:season 2 0.1138033 0.03106555 1.306387 0.230
Residual 69 3.0053980 0.82040127 NA NA
Total 74 3.6633269 1.00000000 NA NA
pairwise.adonis2(beta_q1p$S ~ type*season, data=sample_metadata_row, perm = 999)
$parent_call
[1] "beta_q1p$S ~ type * season , strata = Null , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2      F Pr(>F)    
type         1  0.29577 0.10235 6.3336  0.001 ***
season       1  0.19224 0.06652 4.1165  0.004 ** 
type:season  1  0.06695 0.02317 1.4337  0.205    
Residual    50  2.33495 0.80797                  
Total       53  2.88992 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)  
type         1  0.12582 0.04567 2.3870  0.048 *
season       1  0.09019 0.03274 1.7110  0.126  
type:season  1  0.06174 0.02241 1.1713  0.289  
Residual    47  2.47750 0.89919                
Total       50  2.75526 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df SumOfSqs      R2      F Pr(>F)    
type         1  0.09601 0.06457 3.2849  0.009 ** 
season       1  0.15237 0.10248 5.2131  0.001 ***
type:season  1  0.04012 0.02698 1.3727  0.248    
Residual    41  1.19834 0.80597                  
Total       44  1.48684 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

Functional

betadisper(beta_q1f$S, sample_metadata$type) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008931 0.0044656 1.4188    999  0.214
Residuals 72 0.226611 0.0031474                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          natural protected urban
natural             0.33800 0.117
protected 0.36301           0.375
urban     0.12598   0.40392      
betadisper(beta_q1f$S, sample_metadata$season) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000085 0.00008522 0.0307    999  0.849
Residuals 73 0.202452 0.00277331                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        autumn spring
autumn          0.849
spring 0.86133       
adonis2(beta_q1f$S ~ type,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_y3a7d7l2xknlttbjczil
term df SumOfSqs R2 statistic p.value
type 2 0.02226324 0.05258339 1.998067 0.196
Residual 72 0.40112594 0.94741661 NA NA
Total 74 0.42338918 1.00000000 NA NA
pairwise.adonis2(beta_q1f$S ~ type*season, data=sample_metadata_row, perm = 999)
$parent_call
[1] "beta_q1f$S ~ type * season , strata = Null , permutations 999"

$natural_vs_protected
            Df SumOfSqs      R2       F Pr(>F)    
type         1  0.01858 0.04818  4.2067  0.103    
season       1  0.13428 0.34814 30.3994  0.001 ***
type:season  1  0.01199 0.03108  2.7135  0.179    
Residual    50  0.22087 0.57261                   
Total       53  0.38572 1.00000                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$natural_vs_urban
            Df SumOfSqs      R2      F Pr(>F)   
type         1 0.014221 0.04503  3.691  0.123   
season       1 0.073060 0.23133 18.962  0.004 **
type:season  1 0.047453 0.15025 12.316  0.012 * 
Residual    47 0.181089 0.57339                 
Total       50 0.315824 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$protected_vs_urban
            Df  SumOfSqs       R2       F Pr(>F)  
type         1 -0.001311 -0.00992 -0.4856  0.745  
season       1  0.014945  0.11305  5.5365  0.089 .
type:season  1  0.007891  0.05969  2.9233  0.167  
Residual    41  0.110674  0.83717                 
Total       44  0.132200  1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

attr(,"class")
[1] "pwadstrata" "list"      

7.2 Beta diversity plots

7.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(season) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = season, fill = season, shape = as.factor(type))) +
    geom_point(size = 4) +
    scale_color_manual(values = season_colors) +
  #scale_shape_manual(values = 1:10) 
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Type of pond")

7.2.2 Neutral

beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
#  filter(!sample %in% c("EHI01629", "EHI01630", "EHI01614", "EHI01662", "EHI01446", "EHI01624")) %>% 
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(type) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = type, shape = as.factor(season))) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values=type_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Season",color='Type of pond')

#### In spring

beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
#  filter(!sample %in% c("EHI01629", "EHI01630", "EHI01614", "EHI01662", "EHI01446", "EHI01624")) %>% 
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
    filter(!region %in% c("Eskoriatza","Villabona")) %>% 
 #   filter(season=="autumn") %>%
      filter(season=="spring") %>%
  group_by(type) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = type)) +#, shape = as.factor(type)
    geom_point(size = 4) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values=type_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) + labs(color='Type of pond') 

7.2.2.1 In autumn

beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
#  filter(!sample %in% c("EHI01629", "EHI01630", "EHI01614", "EHI01662", "EHI01446", "EHI01624")) %>% 
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
    filter(!region %in% c("Eskoriatza","Villabona")) %>% 
   filter(season=="autumn") %>%
  #     filter(season=="spring") %>%
  group_by(type) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = type)) +
  geom_point(size = 4) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values=type_colors)+
  theme_classic() +
  theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Type of pond')

7.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
#  filter(!sample %in% c("EHI01629", "EHI01630", "EHI01614", "EHI01662", "EHI01446", "EHI01624")) %>% 
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(season) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = season, shape = as.factor(type))) +
    geom_point(size = 4) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values=season_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Type of pond",color='Season')
beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
#  filter(!sample %in% c("EHI01629", "EHI01630", "EHI01614", "EHI01662", "EHI01446", "EHI01624")) %>% 
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(region) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = region, shape = as.factor(season))) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values=location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Season",color='Type of pond')+
  geom_text(aes(label = sample), size=4)

7.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
#  filter(!sample %in% c("EHI01629", "EHI01630", "EHI01614", "EHI01662", "EHI01446", "EHI01624")) %>% 
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(type) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = type, shape = season)) +
    geom_point(size = 4) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  scale_color_manual(values=type_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Season",color='Type of pond')